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The fractal dimension of a liquid column is a crucial parameter in several models describing the 
main features of the primary break-up occurring at the interface of a liquid phase surrounded by 
the gas-flow. In this work, the deformation of the liquid phase has been numerically studied. The 
gas-phase is computed as a continuum in an Eulerian frame while the liquid phase is discretized 
in droplets Lagrangian tracked and coupled via the momentum equation with the surrounding gas 
flow. The interface is transported by the flow field generated because of the particle forcing and it 
is numerically computed using the Level-Set method. Finally, the fractal dimension of the interface 
is locally estimated and used as criterion for the model of the primary breakup. 



I. INTRODUCTION 

The break-up of a liquid flow leads the evaporation 
process in a wide range of applications such as in combus- 
tion, exhaust gas post-treatment and other devices. The 
rate of break-up enhance the liquid evaporation since the 
surface/volume ratio increases. The atomization occurs 
according with two main mechanisms: the primary and 
the secondary break-up [l|]. 

Firstly, coherent structures detach from the liquid/gas 
interface up to break in spherical droplets because of ei- 
ther the Rayleigh- Taylor or Kelvin- Hclmholtz instability. 
The instability is locally modulated by waves depending 
on the velocity and the pressure fluctuations induced on 
the smaller length scales because of the turbulence (pri- 
mary break-up). 

Thus, each droplet, occurring in the primary breakup, 
fragments up to evaporate because of the shear stress in- 
duced by the surrounding gas-flow (secondary breakup). 

In the case of low pressure injection, the break-up de- 
pends mainly on the surface energy while the interface 
instability plays a negligible role. The surface energy 
could be estimated using a model based on the fractal 
dimension of the liquid/gas interface (2j. The connection 
between the shape of the liquid column and the atomiza- 
tion process has been pointed out in several experimental 
studies using the fractal dimension 

HI!- The topol- 
ogy of the core region is challenging to investigate using 
experiments since the measurements are affected by the 
resuspended droplets detaching from the liquid-gas inter- 
face. 

Numerical simulation appears to be a promising tool 
for this purpose. In the last two decades there has been 
progress on this topic improving the numerical methods 
for the the solution of the governing equations of each 
phase and for the topological treatment of the fluid-gas 
interface. Following, the the most remarkable techniques 
used in literature have been depicted. 

The Front Tracking Jo], the Volume of Fluid method 
0] and the Level-Set [8[ are the main front-capturing 
methods. 

The Front Tracking method is based on the Lagrangian 



treatment of the interface that is described by ideal par- 
ticle (markers). The markers are connected each other 
and they are moved like an ideal fluid particle. The con- 
nectivity between discrete elements covering the interface 
is not straightforward in three spatial dimensions. The 
drawback of this method depends on the collapsing, the 
stretching and the deformation of some discrete elements. 
These elements should be periodically replaced in order 
to preserve the curvature and the mass of the fluid con- 
fined by the interface. These difficulties vanish if the in- 
terface is described using methods in the Eulerian frame 
as the Volume of Fluid method or the Level-Set method. 

In the Volume of Fluid method each phase is repre- 
sented by the volume fraction of the fluid per each com- 
putational cell and the interface evolves according to the 
advection equation. Despite of the simplicity in describ- 
ing the connectivity, a thickness of the interface is ar- 
tificially introduced. Thus, the curvature and the sur- 
face tension are not well estimated. The advantage of 
this method consists in the accurate mass conservation. 
With the Level-Set method the interface is described by 
the iso-contour of an implicit function </>(x) = 0, defined 
in all the points x of a fixed computational domain Q, 
whose evolution is governed by the advection equation. 
The interface, the outer region and the inner region are 
defined by (/>(x) = 0, 4>(x) > 0, 4>(x) < 0, respectively. 
The implicit function ensure the description of complex 
gas-fluid interface without the stringent conditions com- 
ing from the parametric representation of the surface. 
With the Level-Set method the interface is well defined 
and the surface tension is precisely predicted since de- 
pends on the gradient of the implicit function, only. 

Traditionally, the Level-Set involves the numerical so- 
lution in the Eulerian frame for both the liquid and gas 
phases. This method fails for high Reynolds number be- 
cause the steep gradients of density and viscosity across 
the interface requires a fine computational grid making 
tricky the fast simulation of such flows. Moreover, this 
procedure is prohibitive in the industrial applications, 
such as engine combustion, because the liquid column 
evolve in a wide region far away from the injector. 

The Eulerian-Lagrangian approach is a strategy to 
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make faster the simulations preserving a reasonable ac- 
curacy for the interface tracking. The gas-phase is com- 
puted as a continuum in an Eulerian frame while the liq- 
uid phase is discretized in droplets Lagrangian tracked 
and coupled via the momentum equation with the sur- 
rounding gas flow. The liquid particle has to be smaller 
enough to feel all the velocity fluctuations of the gas- 
phase in order to preserve as much as possible the insta- 
bility effects induced on the liquid phase. 

Once the velocity field has been computed both for the 
gas and liquid phases the interface is transported accord- 
ing to the Level-Set method. The advection velocity of 
the interface is equal to the velocity of the liquid par- 
ticles along the interface. Finally, the fractal dimension 
of the liquid/gas interface is locally estimated using the 
box-counting procedure. 



and the mean curvature is: 

k = (V 2 0)/2. (5) 

The signed distance function turns out to be a good 
choice, since steep and flat gradients as well as rapidly 
changing features are avoided as much as possible. The 
evolution of the interface is governed by the following 
equation: 

& + v-V^ = 0, (6) 

where v is the velocity of the liquid phase computed at 
the interface. Once the interface location is known the 
Navier-Stokes equation are solved in the gas-phase and 
the liquid-phase is Lagrangian tracked. 



II. LEVEL-SET METHOD 



III. NUMERICAL METHOD FOR THE 
GAS-PHASE 



In three space dimension a closed surface separate the 
whole domain f2 into the inside domain Q ~ and the out- 
side domain The border between the inside domain 
and the outside domain is called the interface <9f2. An im- 
plicit interface representation defines the interface as the 
iso-contour of some function </>(x) with x e ft. The in- 
terface, the outer region and the inner region are defined 
by </>(x) = 0, </>(x) > 0, 0(x) < 0, respectively 

The gradient of the implicit function V</> is perpendic- 
ular to the iso-contours of (f>. Therefore, V(p evaluated at 
the interface is a vector that points in the same direction 
as the local unit (outward) normal n to the interface. 
Thus, the unit (outward) normal is: 



V0 

W\ 



(i) 



for points on the interface. 

The mean curvature of the interface k is defined as the 
divergence of the normal k = (V • n)/2 so that k > for 
convex regions, k < for concave regions, and k = for 
a plane. Using the definition of normal vector we obtain: 
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(2) 



Smoothness of the is a dcsidcrablc property espe- 
cially in using numerical approximations. Signed dis- 
tance functions are a subset of the implicit functions 4> 
with the extra condition of V0(x) = 1 enforced: 



</>(x) = mm(x — Xj) x G f2, Xi £ <90 (3) 
Under these assumption the normal to the interface is: 



n = 



(4) 



The integral form of the conservative equations applied 
to the control volume Vq with a boundary surface So and 
the unity normal vector pointing outward to the surface 
So defined by n is: 



v 



dp 
lit 



dV 



pu ■ ndS = 



(7) 



and the momentum equation is: 



dpu 



dV 



puu ■ ndS 
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pndS + / (itV u • ndS 

So J So 

[ P gdV+ [ f Uquid dV (8) 

JVn JVn 



where p, fi are the density and the dynamic viscosity, 
respectively. The flow velocity is u, the gravity vector is 
g , the flow pressure is p, the time is t and the feedback 
forcing induced by the liquid phase is fu qu id- 

The numerical solver (KIVA) of the gas-phase equa- 
tions is based on a finite volume method. Spatial dif- 
ference approximations are constructed by the control- 
volume or integral-balance approach, which largely pre- 
serves the local conservation properties of the differential 
equations. In the finite- volume approximations of KIVA, 
velocities are located at the vertices and the scalar quan- 
tities are located at cell centers. Surface and volume 
integrals are approximated using suitable quadrature for- 
mulae. Volume integrals of gradient or divergence terms 
are converted into surface area integrals using the diver- 
gence theorem. The volume integral of a time derivative 
maybe related to the derivative of the integral by means 
of Reynolds transport theorem. Volume and surface area 
integrals are performed under the assumption that the 
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integrands are uniform within cells or on cell faces. Thus 
area integrals over surfaces of cells become sums over cell 
faces i: 



J FdA = Y^FdA 



(9) 



The effect of turbulence is take into account with a 
standard version of k — e turbulence model modified to 
include liquid-turbulence interaction. Evaporating liq- 
uid phase is represented by a discrete-particle technique 
The momentum exchange is treated by implicit coupling 
procedures to avoid the prohibitively small time steps 
that would otherwise be necessary. Turbulence effects 
on the droplets are accounted. When the time step is 
smaller than the droplet turbulence correlation time, a 
fluctuating component is added to the local mean gas 
velocity when calculating each particles momentum ex- 
change with the gas. When the time step exceeds the 
turbulence correlation time, turbulent changes in droplet 
position and velocity are chosen randomly from analyti- 
cally derived probability distributions for these changes. 
The temporal differencing is performed with a first-order 
scheme. 



IV. NUMERICAL METHOD FOR THE 
LIQUID-PHASE 



of a particle to the surrounding gas velocity field is es- 
timated by the Stokes number St = t/t p , in which r is 
the time scale of the eddy interacting with the particle 
characterized by the relaxation time r p . The time scale 
t = L/U depends on the characteristics length L and 
velocity U of the gas-phase. 

Specifically, the droplet with very small Stokes num- 
ber will simply be a flow tracers St « 1 (lighter par- 
ticle). For increasing the Stokes number, the particle 
trajectory will diverge from the flow stream-line until 
the particle motion will be completely independent from 
the surrounding flow St >> 1 (heavier particle). If the 
Stokes number is order one the particle has the maxi- 
mum interaction with the velocity fluctuation of the gas- 
phase [IH . In Figure [JJ it is qualitatively depicted the 
particle-eddy interaction in which the particle velocity 
approaches quickly to the fluid velocity for decreasing 
the Stokes number. 
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FIG. 1: Qualitatively description of the interaction between 
particle and vortices. St « 1 - lighter, St ~ 0(1) - interme- 
diate, St >> 1 - heavier. 



XXX-stochastic method 

The liquid phase has been modelled like a cloud of 
rigid spherical particles with all the forces applied on 
the particle centroid. The time history of each particle 
is numerically computed (Lagrangian particle tracking) 
using a model for the momentum equation with drag and 
gravity force, only [iljflll]: 



dv 

'~dl 



= lc D ^p\vL-v\(u-v)+m pE (10) 



in which m p is the mass, d p is the diameter and Co is 
the drag coefficient. 

The slight deviation from the Stokes flow has been 
taken into account with the drag correction Cd proposed 
by Shiller and Naumann [l2j for particle with Reynolds 



number Re p < 1000. 



C^|i(l + I^) ; Be, 



pd p \u - v| 



(11) 



We can rearrange the equation (Eq. ITU]) as follows: 



dv 
~dt 



- (l + lRei/ 3 ) (u-v) + 



(i 



P P d 2 p 
18/x 



(12) 



The particle is selectively sensitive to the flow velocity 
fluctuations induced by the turbulence. This response 



Particle behavior depending on the Stokes number 
can analytically depicted in the limit of one-dimensional 
Stokes flow (Re p << 1), neglecting the gravity for a 
particle dispersed in a gas-phase with constant velocity 
u = uq. Under this hypothesis the equation 1 121 in dimen- 
sionless form read: 



where 



dv* 
~dt* 
u /U, 



St 



/ * 



(13) 



v/U and t* = tU/L. The 



solution of the previous equation (Eq. [TT| . in the case of 
initial particle velocity equal to zero, has been plotted in 

Fig.LH 



1 



-t* /St 



(14) 



In this scenario, the particle diameter, in the La- 
grangian model of the liquid phase, should be small 
enough in order to make each particle sensitive to the 
velocity fluctuation of the smallest vortical scale. In the 
turbulent regime the smallest vortical scale is defined by 
the dissipative length scale rj, the velocity Uk and the 
time Tfc . Following the Kolmogorov theory the large scale 
(U, L, t) and the dissipative scale (Uk, Lk, Tk) are con- 
nected according to the following relations: 

^(^) 1/4 , U k = (ue)W r fe ^) 1/2 ; (15) 
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FIG. 2: Parametric solution of Eq. [T3]for Uq — 1 and three 
Stokes numbers: St = 0.05, St = 0.1, St = 0.2 



in which e = U 3 /L assess the energy in the whole do- 
main. Thus, the equation 1151 read: 

l = Re-V\ ^ = Re~ 1 / 4 T ± =Re -m. ( 16 ) 
L U t 

where Re = (UL)/v is the Reynolds number evaluated 
at the large scale. 

The Stokes number based on the Kolmogorov scales 
Stk is: 

Stk = 2k = 2kZ. = stRe 1 ' 2 (17) 

Tk T T k 

Thus, if Stk << 1 each particle of the liquid phase 
takes into account the instability effects occurring at the 
liquid/gas interface since the liquid particle is sensitive to 
the velocity fluctuation in the whole range of the turbu- 
lent scales occurring in the computational domain. For 
this reason, the particle diameter should be satisfy the 
equation [TH1 

St k « 1 => StRe 1 ' 2 « 1 =S> 



=> <L <<3\/2- L- [J^-Re- 3 /* (18) 
V Pp 

V. RESULTS 

The time and the spatial evolution of the gas/liquid 
interface has been estimated using the Level-Set method. 
According to the Level-Set method, the evolution of the 
interface is governed by the advection equation [SJ The 



initial condition of Eq. [6] is R s = di/2, in which R a 
is the radius of curvature of the liquid phase and rf,/2 
is the radius of the injector. The liquid particles are 
injected into the gas phase within a prescribed maximum 
angle on. In order to enforce both the angle and the 
injector diameter di the particles have been dispersed at 
the following distance from the orifice: 

*-!•(!) a.) 

The center of the sphere, with radius R s , is defined 
in a Cartesian frame by (x s , y s , z s ), where y s = yt and 
x s = Xi, while z s depends on R s and on as follows: 




FIG. 3: Initial level-set iso-contour. 

The numerical simulations have been carried out using 
an injector with diameter di = 1.58 • 10 -4 \m] crossed by 
a liquid with velocity U — 27 [m/s] and density p p = 10 3 
[Kg/m 3 ]. For such flow the Reynolds number is equal to 
Re — U ■ di/n = 4.310 3 . According to equation ITS1 the 
particle diameter should be d p << 6.9 • 10 -9 [m], thus 
the liquid jet has been discretized using droplets with 
diameter d p = 10 -9 [m]. 

The liquid jet is organized in an elongated structure 
with inhomogeneous shape modulated by the velocity 
fluctuations leading the instability at the liquid/gas in- 
terface. The degree of complexity of the interface has 
been measured using the fractal dimension. The fractal 
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FIG. 4: Interface of the liquid jet. 



dimension is the control parameter in several breakup 
models since it is the integral measure of the distortion 
of the liquid jet 0]. The fractal dimension can be re- 
garded as the general case of the Eulerian dimension. 
If we take an object residing in Euclidean dimension D 
and reduce its linear size by 1/r in each spatial direction, 
its measure (length, area, or volume) would increase to 
N = r D times the original. This is pictured in the next 
figure. 



D=2 
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FIG. 5: Examples of figures with integer dimen- 
sion. Lecture notes of Fluid Dynamics by R. Verzicco 
|http: / /imel 15.poliba.it /verzicco/] 



We consider N — r , take the log of both sides, and 



get 

log(N) = Dlog(r) (22) 

If we solve for D. 

D = log(N)/log(r) (23) 
In the Eulerian frame D is an integer number. Fractals, 
which are irregular geometric objects, have a fractional 
dimension. This theoretical criterion is not straightfor- 
ward to compute the fractal dimension in numerical com- 
putation. For this reason, the computation of the fractal 
dimension of the liquid jet is based on the box-counting 
procedure. To calculate the fractal dimension using the 
box-counting algorithm the liquid/gas interface has been 
placed on a grid. Then, the number of blocks crossed by 
the surface has been stored and iteratively computed for 
resizing grids. The fractal dimension D is the slop of the 
best fit in a log-log scale plane of the number of boxes 
N(s) as function of the size s of the resized grids. The 
box-counting method is widely used since it can measure 
geometry that are not self-similar. 
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FIG. 6: Best fit of the fractal dimension using the box count- 
ing method. 



The slope of the best fit is equal to D = 2.605 and 
it represent the degree of complexity of the liquid jet 
surface. 



VI. SUMMARY AND CONCLUSIONS 
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